%  Written by Walid Boumerdassi and Etienne Collange
%  ECE Georgia Tech
%  Team Space Busters

clear all;
EXIT=0; % Compute the simulated EXIT chart only if EXIT=1
        % If EXIT=0, it means that we are doing a BER analysis on range of
        % SNR. 

BER_IT=0; % When BER_IT=1, it enables a plot of the BER 
         % with respect to the number of iterations
        
nb_iter=6;  % The number of iteration, it need to be adjusted 
            % depending on the SNR because if we take a number of iteration 
            % too high with a high SNR then no error occurs and we stay
            % stuck in the while
            
if EXIT==1 | BER_IT==1
     EbN0_db_vect=0.2;% We compute the trajectory  and/or the Ber_it
                      % only for one SNR
else
     EbN0_db_vect= [0.5:0.1:1];% Different values of SNR that will be used to evaluate the evolution of the BER
    
end
            
% generator polynomials
g=[1 1 1;1 0 1];
% generator polynomials parameters
[n,K] = size(g);
% number of states
m = K-1;
%Code Rate
R=1/K;

mes_size=10000;% Number of random bits generated by block
block_size = mes_size+m;%We need to had m bits to terminate the trellis

%Determine the interleaver by the vector of new positions : pos
[random, pos] = sort(rand(1,block_size));

i_ber=0; %Indice for the BER

ber=zeros(1,length(EbN0_db_vect));
berIT=zeros(length(EbN0_db_vect),nb_iter);
nb_errIT=zeros(1,nb_iter);

for i_ber=1:length(EbN0_db_vect)
    EbN0_db=EbN0_db_vect(i_ber);

    nb_err= 0;

    if (EbN0_db<0.7)
        nb_err_lim= 100/(EbN0_db);
    else
        nb_err_lim=25;
    end
%     if EbN0_db==0
%         nb_err_lim= 100;
%            
%     else
%         if EXIT==1 % When the code is run to plot the trajectory of the decoder
%                    % then the number of errors required does not need to be
%                    % high because the BER is not computed
%            %nb_err_lim=10/EbN0_db;
%            nb_err_lim=2;
%              
%         else
%             nb_err_lim=2;
%            %nb_err_lim=100/(EbN0_db); % The number of errors required
%                                      % to compute the BER is adjusted with
%                                      % the SNR in order to save time in
%                                      % the simulation
%         end
%     end
    nb_mes=0;

    IA1mean=zeros(1,nb_iter);
    IE1mean=zeros(1,nb_iter);
    IA2mean=zeros(1,nb_iter);
    IE2mean=zeros(1,nb_iter);
    
    while nb_err<nb_err_lim % To compute the BER on a significant number of blocks
        % encoder input source code
        InMess = round(rand(1, mes_size));
        %InMess=[1 0 0 0 1 1];
        %InMess=zeros(1,1000);

        %----------------------- Encoder ----------------------------------------


        TurboEncoderOutputs = TurboEncode(g, InMess, pos);

        %interleave x, p1, and p2 in one vector: [x_0 p1_0 p2_0 x_1 p1_1 p2_1 ...]
        encodedBits= zeros(1,3*length(TurboEncoderOutputs(1,:)));
        for i=1:block_size
            encodedBits(3*i-2:3*i) = (TurboEncoderOutputs(:,i)');
        end

        % Transform the bits in symboles: 0->-1 and 1->+1
        encodedSymb = 2*encodedBits - 1;

        %----------------------- Channel ----------------------------------------

        EbN0=10^(EbN0_db/10);

        sigmaB_2= 1/(2*R*EbN0);

        noise= randn(1,length(encodedSymb)).*sqrt(sigmaB_2);

        noisySymb= encodedSymb + noise; %Add the noise to the encoded symboles

        % Turbo decoder
        %-----------------------------------------------------
        Lc=2/sigmaB_2;

        % Reorder x, p1 and p2 in a matrix of three rows
        noisyDecIn = [ noisySymb(1:3:end);
            noisySymb(2:3:end);
            noisySymb(3:3:end)];

        extrinsic1 = zeros(block_size,1);
        extrinsic2 = zeros(block_size,1);
        apriori1 = zeros(block_size,1);
        apriori2 = zeros(block_size,1);

        xp= noisySymb(1:3:end);

        InDec1=zeros(1,2*block_size);
        InDec1(1:2:end)= xp;
        InDec1(2:2:end)= noisySymb(2:3:end);


        InDec2=zeros(1,2*block_size);
        InDec2(1:2:end)= xp(pos);
        InDec2(2:2:end)= noisySymb(3:3:end);
        Le1=zeros(block_size,nb_iter);
        Le2=zeros(block_size,nb_iter);
        IA1=zeros(1,nb_iter);
        IE1=zeros(1,nb_iter);
        IA2=zeros(1,nb_iter);
        IE2=zeros(1,nb_iter);
        for j=1:nb_iter

            % First decoder
            apriori1(pos) = extrinsic2;

            %-----------Computation of mutual information between apriori1 and X
            if EXIT==1
                LaPlus= apriori1(find(TurboEncoderOutputs(1,:)>0.5));
                LaMinus= apriori1(find(TurboEncoderOutputs(1,:)<0.5));

                pAp=hist(LaPlus,min(apriori1):1:max(apriori1));
                pAp=pAp/sum(pAp);% pAp(x) = pA(x|X=1)

                pAm=hist(LaMinus,min(apriori1):1:max(apriori1));
                pAm=pAm/sum(pAm);% pAm(x) = pA(x|X=-1)

                IA1(j)=1/2*(sum(pAp(find(pAp>0)).*log2(2*pAp(find(pAp>0))./...
                    (pAp(find(pAp>0))+pAm(find(pAp>0)))))...
                    +sum(pAm(find(pAm>0)).*log2(2*pAm(find(pAm>0))./...
                    (pAp(find(pAm>0))+pAm(find(pAm>0))))));
            end
            %-----------Output of First decoder
            [LMAP1 extrinsic1] = TurboDecode1(g, noisyDecIn,Lc, block_size, apriori1);

            %-----------Computation of mutual information between extrinsic1 and X
            if EXIT==1
                extrPlus= extrinsic1(find(TurboEncoderOutputs(1,:)>0.5));
                extrMinus= extrinsic1(find(TurboEncoderOutputs(1,:)<0.5));

                pEp=hist(extrPlus,min(extrinsic1):1:max(extrinsic1));
                pEp=pEp/sum(pEp);% pEp(x) = pE(x|X=1)

                pEm=hist(extrMinus,min(extrinsic1):1:max(extrinsic1));
                pEm=pEm/sum(pEm);% pEm(x) = pE(x|X=-1)

                IE1(j)=1/2*(sum(pEp(find(pEp>0)).*log2(2*pEp(find(pEp>0))./...
                    (pEp(find(pEp>0))+pEm(find(pEp>0)))))...
                    +sum(pEm(find(pEm>0)).*log2(2*pEm(find(pEm>0))./...
                    (pEp(find(pEm>0))+pEm(find(pEm>0))))));
            end
            %-----------Stocking of extrinsic1 in Le1
            Le1(:,j)=extrinsic1;

            %Second decoder
            apriori2 = extrinsic1(pos);

            %-----------Computation of mutual information between apriori2 and X
            if EXIT==1
                LaPlus= apriori2(find(TurboEncoderOutputs(1,pos)>0.5));
                LaMinus= apriori2(find(TurboEncoderOutputs(1,pos)<0.5));

                pAp=hist(LaPlus,min(apriori2):1:max(apriori2));
                pAp=pAp/sum(pAp);% pAp(x) = pA(x|X=1)

                pAm=hist(LaMinus,min(apriori2):1:max(apriori2));
                pAm=pAm/sum(pAm);% pAm(x) = pA(x|X=-1)

                IA2(j)=1/2*(sum(pAp(find(pAp>0)).*log2(2*pAp(find(pAp>0))./...
                    (pAp(find(pAp>0))+pAm(find(pAp>0)))))...
                    +sum(pAm(find(pAm>0)).*log2(2*pAm(find(pAm>0))./...
                    (pAp(find(pAm>0))+pAm(find(pAm>0))))));
            end
            %-----------Output of Second decoder
            [LMAP2 extrinsic2] = TurboDecode2(g, noisyDecIn,Lc, block_size, apriori2, pos);

            %-----------Computation of mutual information between extrinsic2 and X
            if EXIT==1
                extrPlus= extrinsic2(find(TurboEncoderOutputs(1,pos)>0.5));
                extrMinus= extrinsic2(find(TurboEncoderOutputs(1,pos)<0.5));

                pEp=hist(extrPlus,min(extrinsic2):1:max(extrinsic2));
                pEp=pEp/sum(pEp);% pEp(x) = pE(x|X=1)

                pEm=hist(extrMinus,min(extrinsic2):1:max(extrinsic2));
                pEm=pEm/sum(pEm);% pEm(x) = pE(x|X=-1)

                IE2(j)=1/2*(sum(pEp(find(pEp>0)).*log2(2*pEp(find(pEp>0))./...
                    (pEp(find(pEp>0))+pEm(find(pEp>0)))))...
                    +sum(pEm(find(pEm>0)).*log2(2*pEm(find(pEm>0))./...
                    (pEp(find(pEm>0))+pEm(find(pEm>0))))));
            end
            %-----------Stocking of extrinsic2 in Le2
            Le2(:,j)=extrinsic2;
            
            if BER_IT==1
                XhatIT=(LMAP1(1:end-2)>0)';
                diffIT=InMess-XhatIT;
                nb_errIT(j)=nb_errIT(j)+sum(abs(diffIT));
            end

        end %for

        %compute the average IAs and IEs
        IA1mean=IA1mean+IA1;
        IE1mean=IE1mean+IE1;
        IA2mean=IA2mean+IA2;
        IE2mean=IE2mean+IE2;


        Xhat=(LMAP1(1:end-2)>0)';
        diff=InMess-Xhat;
        nb_err=nb_err+sum(abs(diff))
        nb_mes=nb_mes+1;
    end%while
    
    ber(i_ber)=nb_err/mes_size/nb_mes
    
    if BER_IT==1
        berIT(i_ber,:)=nb_errIT/mes_size/nb_mes
    end
end% for EbN0-db

if EXIT==1
    IA1mean=IA1mean/nb_mes;
    IE1mean=IE1mean/nb_mes;
    IA2mean=IA2mean/nb_mes;
    IE2mean=IE2mean/nb_mes;
    j=0;
    for i=1:length(IA1mean)
        j=j+1;
        traj(j)=IE1mean(i);
        indices(j)=IA1mean(i);
        j=j+1;
        traj(j)=IA2mean(i);
        indices(j)=IE2mean(i);
    end
    plot(indices,traj,'k.-')
elseif BER_IT==1
    figure;
    plot(berIT(1,:));
    xlabel('Number of iterations');
    ylabel('BER');
else
    figure;
    semilogy(EbN0_db_vect,ber,'r*-')
    xlabel('EbN0_db')

    ylabel('ber')
    hold on 
    BPSK_WO=1/2.*erfc(sqrt(10.^(EbN0_db_vect./10)));
    semilogy(EbN0_db_vect, BPSK_WO,'b*-')
    hold off
    
end